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1.  SUMMARY 


We  are  investigating  generation  of  complex  seismic  waves  by  explosions  in  media  with  realistic 
3D  heterogeneity,  including  topography  and  tectonic  strain  release.  We  have  developed  a  3D 
nonlinear  finite  element  code  designed  for  calculation  of  explosions  in  3D  heterogeneous  media 
with  well-tested  material  models  and  including  gravity;  and  implemented  interface  codes  using 
the  representation  theorem  that  propagate  the  numerical  solution  from  the  source  region  to 
regional  and  teleseismic  distances. 

We  have  run  calculations  of  the  NPE  with  and  without  topography,  and  extended  the  calculation 
to  400  km  using  the  representation  theorem.  The  calculation  shows  that  the  effect  of  topography 
can  be  modeled  and  produces  effects  consistent  with  the  data.  Near  field  tangential  motion  is 
caused  primarily  by  topographic  gradients.  The  effect  of  near-source  topography  on  subsurface 
ground  motion  is  primarily  due  to  off-axis  propagation.  The  effects  of  near-source  topography  on 
far-field  seismic  radiation  are  small. 

We  perform  calculations  of  the  nuclear  test  SHOAL  with  and  without  tectonic  strain  release.  The 
calculation  shows  that  tectonic  strain  release  consistent  with  the  local  stress  state  can  generate 
near  field  and  regional  signals,  including  long  period  surface  waves,  similar  to  those  observed. 
We  find  that  tectonic  release  causes  small  changes  to  the  far- field  P-wave  waveform,  but  has 
very  little  effect  on  far-field  P-wave  amplitudes.  It  does  generate  SH  waves  not  present  in  the 
calculation  without  tectonic  release.  This  implies  that  while  tectonic  release  may  significantly 
affect  Ms,  the  effect  on  mb  will  be  small.  We  also  find  that  nonlinear  deformation  caused  by  free 
surface  interaction  increases  the  Mzz  component  of  the  source,  adding  a  CLVD  component. 
However,  the  free  surface  interaction  also  increases  the  Mxx  and  Myy  components  and  the  net 
effect  on  surface  wave  amplitudes  is  small. 

We  use  available  information  on  structural  heterogeneity  to  develop  3D  models  that  have  the 
statistical  characteristics  of  the  upper  few  kilometers  of  the  earth’s  crust.  We  modified  the  Shoal 
earth  model  to  have  these  statistical  characteristics.  We  ran  a  series  of  calculations  first  varying 
the  elastic  properties  of  the  structure,  but  keeping  strength  the  same,  and  second  modifying  the 
strength  to  vary  as  a  function  of  shear  modulus.  The  extra  degree  of  freedom  provided  in  three 
dimensions  makes  generation  of  shear  waves  easier  than  with  the  restriction  to  axisymmetric 
geometry.  All  of  calculations  with  the  heterogeneous  structure  generated  shear  waves  not  present 
in  the  initial  plane-layered  structure,  however  the  effect  is  much  stronger  with  strength  variation. 

2.  INTRODUCTION 

This  project  follows  a  previous  project  on  3D  numerical  modeling  and  wave  propagation.  In  that 
project,  we  did  the  following: 

1 .  Developed  a  3D  nonlinear  finite  element  code  designed  for  calculation  of  explosions  in 

3D  heterogeneous  media  with  well-tested  material  models  and  including  gravity; 

2.  Implemented  interface  codes  using  the  representation  theorem  that  propagate  the 

numerical  solution  from  the  source  region  to  regional  and  teleseismic  distances; 

3.  Used  the  code  to  evaluate  the  effect  of  near  source  heterogeneity  including  non-linear 

material  property  variations,  elastic  variability  and  topography. 
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The  reason  that  3D  calculations  are  important  for  understanding  shear  wave  generation  is  that 
symmetry  constraints  imposed  by  ID  and  2D  calculations  act  to  suppress  shear  waves. 
Imposition  of  2D  axisymmetry  requires  that  the  motion  generated  by  the  explosion  be  identical 
in  all  horizontal  directions  through  propagation  of  the  shock  wave,  nonlinear  deformation  of  the 
surrounding  material,  and  rebound  to  a  final  state.  Since  the  motion  is  never  identical  in  all 
directions,  shear  waves  will  always  be  generated,  and  in  fact  shear  waves  are  observed  from  all 
underground  explosions.  In  this  project  we  are  attempting  to  develop  realistic  representative 
models  for  3D  near-source  heterogeneity,  followed  by  3D  explosion  calculations  in  these 
models.  We  have  performed  calculations  for  two  specific  events:  the  Non-Proliferation 
Experiment  (NPE)  and  the  Shoal  underground  nuclear  test.  For  the  NPE,  we  included  the 
topography  above  and  near  the  explosion.  For  Shoal,  we  included  tectonic  prestress 
corresponding  to  the  stress  state  for  that  location  in  Nevada. 

3.  TECHNICAL  APPROACH 
3.1  The  CRAM3D  Code 

CRAM  3D  is  an  explicit  three-dimensional  Lagrangian  finite  element  code  designed  to  run  on 
multiple  processors  (Stevens  and  Xu,  2010;  Stevens  et  al,  2011;  Stevens  et  al,  2014).  For  an 
explosion  simulation,  the  cavity  is  placed  near  the  center  of  the  grid  and  is  enclosed  by  a  spider 
grid  which  facilitates  applying  the  pressure  boundary  condition  and  rezoning  elements  (Figure 
1).  The  well-tested  nonlinear  material  models  from  CRAM  2D  have  been  implemented  in 
CRAM  3D.  The  code  includes  gravity  and  so  includes  the  important  effects  that  result  from 
variation  of  overburden  pressure  with  depth.  Gravitational  equilibrium  is  established  by  running 
an  initial  calculation  with  no  source,  followed  by  a  second  calculation  including  the  explosion 
source.  We  have  also  implemented  the  capability  to  include  tectonic  prestress  in  the  calculations. 


Figure  1.  The  CRAM  3D  finite  element  outer  grid  (left)  is  rectangular.  The  inner  grid  (center)  is 
shaped  to  match  the  shape  of  the  explosion  shock  wave.  CRAM2D  uses  a  similar  axisymmetric 
spider  grid  (right)  in  the  region  around  the  explosion. 
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3.2  Propagation  with  the  Elastodynamic  Representation  Theorem 

The  representation  theorem  allows  us  to  perform  arbitrarily  eomplex  nonlinear  ealeulations  in  the 
source  region,  and  then  propagate  them  with  an  appropriate  Green’s  function.  The  representation 
theorem  is  exact.  That  is,  no  matter  how  complex  the  3D  motion  is  on  the  source  region 
boundary,  it  will  be  correctly  propagated  by  the  representation  theorem.  The  only  exception  is 
that  it  will  not  calculate  the  interaction  of  backscattered  waves  reflected  from  outside  the  source 
region  with  complexities  of  the  source  region. 

In  the  three-dimensional  numerical  finite  difference  calculations,  we  save  displacements  and 
stresses  (or  velocities  and  stress  rates)  due  to  the  seismic  source  on  a  monitoring  surface  on  the 
boundary  of  a  rectangle  (5  planar  surfaces,  excluding  the  upper  surface),  and  calculate  Green’s 
functions  from  each  point  on  the  monitoring  surface  to  the  receiver  and  so  the  synthetic 
seismogram  at  the  receiver  point  X  outside  of  the  monitoring  surface  is  obtained  by  integrating 
over  the  monitoring  surface 

Ui(X)  =  I  {g5(^;  X)  .  rf  (0  -  uf  (0  HC  X)7i,]dA  (1) 

i>M 

in  the  frequency  domain,  where  Gj{^;X)  and  Sj,^(^;X)  are  the  Green’s  function  and  the  stress 

tensor  on  the  monitoring  surface  due  to  a  unit  impulsive  force  at  X  in  direction  i,  is  the 

traction  on  the  monitoring  surface  due  to  the  seismic  source,  u  is  the  displacement  on  the 
monitoring  surface,  and  n  is  the  normal  to  the  monitoring  surface.  The  operator  *  denotes 
convolution  and  the  summation  convention  is  assumed. 

We  use  a  plane-layered  Green’s  function  outside  the  source  region.  The  Green’s  functions  for  the 
complete  seismograms  are  derived  from  an  algorithm  based  on  the  work  of  Luco  and  Apsel 
(1983)  and  Apsel  and  Luco  (1983).  The  technique  used  for  surface  waves  is  similar  to  the 
method  of  Bache  et  al.  (1982).  The  Green’s  functions  for  body  waves  are  generated  by  a 
procedure  similar  to  that  described  by  Bache  and  Harkrider  (1976)  using  a  saddle  point 
approximation  to  calculate  a  far-field  plane  wave  for  a  given  takeoff  angle  from  a  source  in  a 
plane-layered  medium.  Although  the  full  waveform  Green’s  functions  generate  the  complete 
waveform,  the  other  Green’s  functions  provide  additional  insight  into  the  source  and  waveform 
generation. 

3.3  Incorporation  of  Tectonic  Prestress 

The  effect  of  tectonic  prestress  on  explosion-generated  surface  waves  has  been  discussed  in 
many  prior  studies  (e.g.,  Archambeau,  1972;  Rygg,  1979;  Toksoz  and  Kehrer,  1972).  In  most  of 
these  studies  tectonic  release  was  modeled  as  superposition  of  a  tectonic  source  described  by  a 
double  couple,  multipole  or  moment  tensor  source,  plus  a  point  explosion  source.  The  size  of  the 
tectonic  source  was  determined  by  comparison  with  the  observed  Love  waves  and  Rayleigh 
wave  radiation  pattern.  Day  et  al.  (1987)  first  attempted  to  perform  numerical  modeling  of 
tectonic  release  through  an  axisymmetric  calculation  of  the  explosion  Piledriver.  Murphy  et  al. 
(2011,  2013)  performed  similar  calculations  for  the  North  Korean  explosions.  To  the  best  of  our 
knowledge  no  one  has  previously  performed  numerical  calculations  for  a  three-dimensional 
stress  field. 
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The  first  step  in  performing  a  calculation  of 
tectonic  release  is  to  determine  the  stress 
state  of  the  earth.  While  the  stress  state  varies 
between  regions,  considerable  work  has  also 
been  done  on  this  subject  (e.g.,  Zoback, 

2010).  In  general  the  vertical  stress  is  equal 
to  the  overburden  weight  of  the  material 
above  at  any  given  point.  The  horizontal  Figure  2.  Stress  regimes  vary  by  region, 

stresses  may  be  larger  or  smaller  than  this  Horizontal  stresses  may  be  substantially  larger 

value  up  to  the  point  where  failure  due  to  qj-  smaller  than  the  vertical  stress, 
frictional  sliding  or  rock  fracturing  relieves 

the  stress.  In  most  of  the  earth’s  crust  tectonic  stresses  are  applied  continuously,  and  so  the 
prestress  is  in  general  close  to  equilibrium  with  frictional  sliding.  The  ratio  of  maximum  to  the 
minimum  principal  stress,  both  reduced  by  pore  pressure,  is  given  by 


Sv=OvBr 

burden 

Weight 


Normal  5hmin<5v 
y  5Hmax<5v 


Strike  Slip  Shmin<Sv 
SHmax>Sv 


Reverse  Shmin>Sv 
SHmax>Sv 


(2) 


where  p  is  the  coefficient  of  friction,  Pp  is  the  pore  pressure,  Smax  is  the  maximum  of  the  vertical 
stress  Sv  and  the  maximum  horizontal  stress  SHmax  and  Smin  is  the  minimum  of  Sy  and  the 
minimum  horizontal  stress  Shmin.  For  a  typical  coefficient  of  friction  of  0.6,  this  ratio  is  3.1,  so 
the  difference  between  principal  stresses  can  be  substantial  and  the  effect  on  seismic  waves 
potentially  large.  In  our  calculations,  we  calculate  the  overburden  weight  to  determine  t'he 
vertical  stress,  and  then  modify  the  horizontal  stresses  to  some  fraction  of  the  frictional  limit. 
This  is  the  initial  stable  state  of  the  calculation  prior  to  introduction  of  the  explosion.  Note  that 
although  the  vertical  stress  is  still  equivalent  to  the  overburden  weight,  the  pressure  is  not,  and  it 
may  be  either  increased  or  reduced  by  the  tectonic  stresses.  Since  material  strength  increases 
with  pressure,  this  also  can  substantially  affect  the  seismic  source.  In  general,  normal  faulting 
regimes  will  amplify  seismic  signals,  while  reverse  faulting  regimes  will  decrease  seismic 
signals;  strike-slip  regimes  could  do  either  (Figure  2). 


3.4  Incorporation  of  Topography 

CRAM3D  incorporates  topography  by  gradually  distributing  offsets  of  the  grid  from  top  to 
bottom  so  that  change  in  cell  dimensions  is  gradual  rather  than  abrupt  at  the  surface.  The  main 
difficulty  in  implementing  topography  is  establishing  gravitational  equilibrium.  The  variation  in 
overburden  pressure  with  depth  is  very  important  in  nonlinear  calculations  because  strength  is 
pressure  dependent  and  so  material  becomes  much  weaker  at  shallow  depths.  Because  of  this  we 
need  to  bring  the  grid  into  gravitational  equilibrium  before  running  an  explosion  calculation.  As 
mentioned  above,  this  is  done  by  starting  with  an  approximate  solution  for  overburden  pressure 
in  the  grid  and  then  running  an  initial  calculation  without  the  explosion.  While  this  can  be 
accomplished  relatively  easily  for  a  plane-layered  structure,  bringing  the  grid  to  gravitational 
equilibrium  in  a  structure  with  significant  topography  has  proven  to  be  difficult.  While  the 
calculation  will  eventually  come  to  equilibrium,  adding  topography  is  the  equivalent  of  dropping 
a  mountain  onto  the  free  surface,  and  so  it  takes  quite  a  long  time  to  come  to  equilibrium.  We 
have  experimented  with  several  approximate  solutions  (e.g.,  Liu  and  Zoback,  1992),  but  have  not 
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as  yet  found  a  solution  that  allows  the  equilibrium  run  to  be  eompleted  without  taking 
considerable  time.  The  calculated  displacements  caused  by  topography  close  to  the  NPE 
calculation  are  shown  in  Figure  3. 


r^nrtti  Range  (m)  Nbnh  Range  (m)  Nbnh  Range  (m) 


Figure  3.  From  left  to  right,  vertical,  eastward  and  northward  displacements  caused  by 
topography  along  a  slice  204  meters  from  the  NPE  shot  point. 

4.  RESULTS  AND  DISCUSSION 

4.1  Calculation  of  the  NPE 

As  our  first  calculation  of  a  historical  event  with  3D  heterogeneity,  we  model  the  Non- 
Proliferation  Experiment  (NPE),  including  the  surface  topography.  The  NPE  was  a  chemical 
explosion  with  yield  equivalent  to  one  kiloton  of  TNT.  The  explosive  is  1.315  kilotons  of  a  50/50 
ANFO/emulsion  mix  in  a  cylindrical  cavity  centered  at  389  m  depth,  7.7  m  in  radius 
(horizontally)  and  5.2m  in  height  (vertically).  We  previously  modeled  this  event  in  a  two- 
dimensional  axisymmetric  calculation  (Stevens  et  al,  2004),  and  initially  used  the  same  material 
properties,  which  were  derived  from  Rimer  et  al  (1994)  for  the  3D  calculation.  The  material 
geology  at  the  NPE  site  is  based  on  that  of  the  nearby  Misty  Echo  event  and  consists  of  four 
layers.  Layers  1  and  4  are  non-porous,  and  layers  2  and  3  have  porosities  of  3%  and  0.5% 
respectively.  Comparison  of  the  initial  calculations  with  recordings  from  this  event  showed  that 
calculated  arrival  times  were  early  relative  to  the  data,  so  porous  moduli  in  the  material  were 
reduced  so  that  P-velocities  were  more  consistent  with  the  observed  arrival  times,  similar  to 
material  velocities  used  by  Kamm  and  Bos  (1995). 
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Figure  4.  Topography  around  the  area  of  the  NPE  calculation  (top  left)  and  the  topography  on  the 
CRAM  grid  surface  (top  right).  Lower  figure  shows  a  Google  Earth  image  with  a  line  due  north 
from  directly  above  the  explosion  to  the  edge  of  the  grid. 


Figure  4  shows  the  actual  topography  near  the  NPE  (National  Elevation  Dataset  1/3  Arc-Second 
http://seamless.usgs.gov/nedl3.php).  For  this  calculation,  we  modeled  the  area  within  650  meters 
of  the  explosion.  Even  that  relatively  small  area  has  substantial  topography,  with  elevation 
differences  of  up  to  140  meters.  The  strongest  topographic  changes  are  along  a  canyon  to  the 
north  and  a  ridge  to  the  south.  The  calculation  area  is  tapered  at  the  edges  to  bring  it  to  a  constant 
elevation  at  and  outside  the  monitoring  surface  that  is  used  to  propagate  the  calculation  as 
described  above. 

Figure  5  (left)  shows  the  pressure  field  output  from  the  calculation  with  topography  at  several 
times.  The  upgoing  shock  wave  interacts  with  the  surface  in  a  complicated  way  causing 
distortion  of  the  reflected  waves.  Figure  5  (right)  shows  the  region  of  nonlinear  deformation  and 
tensile  cracking.  Nonlinear  deformation  near  the  explosion  is  approximately  spherical,  but 
elongated  and  slightly  offset  vertically.  The  region  of  cracking  is  confined  to  near  the  free 
surface.  This  explosion  is  deeply  overburied,  so  there  is  less  asymmetry  than  in  a  normally 
buried  explosion. 

Figure  6  shows  the  calculated  and  observed  velocity  waveforms  at  three  points  on  the  free 
surface.  Spall  is  indicated  by  a  linear  trending  velocity  with  a  slope  of  g,  and  is  readily  apparent 


Approved  for  public  release;  distribution  is  unlimited. 
6 


in  the  ground  zero  waveform  and  also  visible  in  the  more  distant  surface  stations.  The  match 
between  calculated  and  observed  waveforms  is  quite  good. 


Figure  5.  Left:  Pressure  field  along  a  vertical  slice  120  meters  from  the  origin  at  6  times:  0.05, 
0.1,  0.2,  0.3,  0.5  and  1.0  seconds.  Note  the  strong  interactions  between  the  shock  wave  and  both 
the  free  surface  topographic  features  and  material  layer  boundaries.  Right:  final  nonlinear 
deformation  along  a  cross  section  through  the  source  showing  regions  of  yielding  and  tensile 
cracking. 


NPE  Near  Surface  Stations 


time  (s) 

Figure  6.  Comparison  of  calculated  and  observed  waveforms  from  the  NPE  on  the  surface  at 
three  locations:  ground  zero,  immediately  above  the  explosion,  and  at  distances  of  264  and  670 
meters.  Note  the  distinct  spall  phase  that  decreases  in  duration  with  distance  in  both  the  data  and 
the  calculation. 
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In  order  to  assess  the  effect  of  topography  distinct  from  other  mechanisms  that  might  affect 
waveforms,  we  also  ran  a  calculation  with  flat  topography.  In  general  the  results  are  not  very 
different.  With  topography  there  is  some  tangential  motion  not  present  in  the  flat  calculation  and 
small  differences  in  waveforms,  but  the  differences  are  not  significant  except  on  the  surface  near 
topographic  gradients.  Figure  7  shows  a  comparison  of  near  field  observed  and  calculated 
waveforms  at  shot  depth.  In  a  uniform  medium  only  the  radial  component  would  be  nonzero. 
However,  the  free  surface  interaction  causes  a  large  arrival  on  the  vertical  component  at  about 
0.7  seconds.  Neither  the  flat  nor  topographic  calculation  predicts  significant  tangential  motion  at 
this  location.  Although  there  is  motion  on  the  tangential  component,  the  clear  P  arrival  and 
strong  resemblance  to  the  vertical  component  suggests  that  the  motion  is  mostly  off-axis  vertical 
or  radial  motion  rather  than  true  tangential  motion. 


T14  -  Okm,  30’  (0  -25  Hz) 


Figure  7.  Shot  level  calculated  and  observed  waveforms  at  a  distance  of  500  meters.  Observed 
tangential  motion  is  most  likely  off-axis  radial  and/or  vertical  motion.  Blue  line:  NPE  data, 
Green  line:  calculated  velocity  without  topography.  Red  line:  calculated  velocity  with 
topography. 

The  strongest  effect  of  topography  is  on  surface  stations,  and  so  we  compared  our  results  with 
the  observed  ground  motion  at  those  locations.  Figure  8  shows  the  calculated  tangential/radial 
component  amplitude  ratio  on  the  surface  above  the  explosion  at  the  locations  where  stations 
were  placed  during  the  event.  Note  that  the  tangential  component  amplitudes  correlate  with 
topographic  gradients.  That  is,  locations  where  the  surface  elevation  changes  rapidly  generates 
tangential  motion.  Figure  9  shows  a  comparison  of  the  observed  and  calculated  tangential/radial 
motion  using  the  first  part  of  the  observed  waveforms  at  each  station.  The  agreement  is  quite 
good,  confirming  that  surface  topography  is  likely  responsible  for  the  tangential  motion 
observed.  Figure  10  shows  a  more  detailed  comparison  of  the  calculated  and  observed  data  over 
a  longer  time  period.  The  agreement  is  reasonably  good.  Spall,  for  example,  is  modeled  quite 
well.  The  main  difference  is  in  the  large  SH  observed  at  later  times  at  some  of  the  stations 
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(bottom  three  stations  in  the  seismogram  plot).  Note  the  strong  resemblance  of  the  phase  to  an 
inverted  direct  P-wave.  There  is  also  a  large  velocity  contrast  in  the  earth  model  just  above  the 
source,  and  there  is  some  slope  to  these  earth  layers  that  we  are  not  modeling  in  our  calculations. 
Our  conclusion  from  this  is  that  the  large  anomalous  “SH”  motion  is  actually  P-wave  that  has 
reflected  off  the  surface,  then  back  up  off  the  reflecting,  dipping  layers  so  that  it  is  no  longer 
vertically  polarized  and  appears  on  the  tangential  component. 
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Figure  8.  Full  topography  (top)  and  tapered  topography  (bottom)  together  with  calculated 
Tangential/Radial  amplitude  ratio.  Tangential  component  data  is  generated  particularly  at 
topographic  gradients. 
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NPE  Stations  on  10m  National  Elevation  Dataset 
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Figure  9.  Tangential/Radial  amplitude  ratio  calculated  with  topography  (top)  and  data  (bottom) 
using  the  early  part  of  the  waveform.  Agreement  is  fairly  good  within  the  region  where 
topography  is  explicitly  modeled  indicating  that  the  early  time  tangential  component  data  on  the 
surface  most  likely  is  caused  by  topography. 
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We  used  the  representation  theorem  to  ealeulate  waveforms  for  eomparison  with  observations  at 
a  number  of  loeations.  We  used  an  earth  model  derived  by  merging  the  local  NPE  model  of 
Rimer  et  al  (1994),  and  the  regional  3D  model  of  Myers  et  al  (2007)  (Table  1).  The  Q  structure 
used  in  the  calculation  is  based  on  Swanger’s  Law,  Q=p/10  where  p  is  shear  velocity.  The 
calculated  and  observed  waveforms  at  ELK  are  shown  in  Figure  11.  The  waveforms  are 
generally  similar,  and  the  relative  P  and  S  amplitudes  are  quite  good.  Arrival  times  are  earlier  in 
the  observations  than  the  calculation,  indicating  that  the  structure  used  is  somewhat  slower  than 
the  actual  earth  structure  between  the  NPE  and  ELK.  Near  field  topography  does  not  have  a 
noticeable  effect  on  the  waveforms  at  this  distance. 


ELK  - 

I 


Figure  11.  Calculated  and  observed  waveforms  from  the  NPE  at  station  ELK  at  a  distance  of  400 
km.  Top:  Observed;  Middle:  calculated  with  topography;  Bottom:  calculated  with  flat 
topography. 


Table  1.  NPE  rej 

?ion  earth  model 

Depth 

Thickness 

Vp 

Vs 

Density 

Q 

(km) 

(km) 

(km/sec) 

(km/sec) 

(gm/cc) 

0.102 

0.102 

3.586 

1.809 

2.200 

180 

0.328 

0.226 

1.824 

1.000 

1.663 

100 

0.704 

0.376 

2.701 

1.357 

1.900 

140 

3.600 

2.896 

3.210 

1.900 

1.680 

190 

3.660 

0.060 

4.500 

2.600 

2.300 

260 

8.100 

4.440 

5.930 

3.500 

2.430 

350 

12.240 

4.140 

5.950 

3.510 

2.453 

351 

14.280 

2.040 

5.970 

3.520 

2.458 

352 

17.280 

3.000 

6.000 

3.540 

2.467 

354 

20.280 

3.000 

6.070 

3.580 

2.487 

358 

23.280 

3.000 

6.200 

3.660 

2.523 

366 

26.220 

2.940 

6.310 

3.720 

2.554 

372 

29.220 

3.000 

6.360 

3.750 

2.568 

375 

34.740 

5.520 

6.410 

3.780 

2.582 

378 

37.260 

2.520 

7.900 

4.400 

3.001 

440 
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4.2  Calculation  of  Shoal 


Shoal  was  a  12.5  kiloton  nuclear  explosion  detonated  near  Fallon,  Nevada  at  367  m  depth 
(Weart,  1965).  There  were  three  near- field  shot  level  recording  stations  located  in  three 
directions  each  at  about  590  meters  from  the  shot  (Figure  12).  The  minimum  principal  stress 
Shmin  in  this  part  of  Nevada  is  in  extension.  It  is  unclear  from  regional  tectonics  whether  the 
maximum  horizontal  stress  SHmax  is  less  than  or  greater  than  the  vertical  stress  Sv.  For  the 
purpose  of  this  calculation,  we  made  SHmax  20%  larger  than  Sv,  and  made  Shmin  95%  of  the 
critical  stress  corresponding  to  frictional  sliding  with  a  coefficient  of  friction  of  0.6.  The  stress 
state  in  the  upper  kilometer  is  shown  in  Figure  13. 


Figure  12.  Location  of  the  Shoal  explosion  and  the  three  near-field  shot  level  recording  stations 
and  the  direction  of  the  local  stress  state.  The  extensional  stress  Shmin  is  the  west-northwest 
direction,  which  is  the  direction  from  the  origin  to  PM-3.  The  second  principal  stress  is  harder  to 
determine  and  in  the  calculation  we  model  it  as  slightly  compressive. 


Stress  State 


Figure  13.  Vertical  and  horizontal  stresses  used  in  the  calculation  shown  vs.  depth.  Dashed  lines 
indicate  the  critical  stress  limited  by  balance  against  friction. 
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Elastic  material  properties  used  in  the  ealculation  were  derived  from  Beers  (1964)  and  are  listed 
in  Table  2.  The  same  material  was  used  throughout  the  grid,  assuming  negligible  porosity  and 
uniform  strength.  Initial  strength  corresponded  to  fraetured  granite,  and  strength  reduction  was 
modeled  using  a  shear  strain  based  shock  damage  model  derived  for  Piledriver  granite  (Stevens 
et  al,  2003).  The  grid  used  in  the  calculation  was  2  km  x  2  km  by  0.8  km  deep  with  10  meter 
spaeing.  The  initial  cavity  was  4.32  meters  in  radius  and  the  inner  grid  spacing  started  at  0.5 
meters,  gradually  increasing  to  the  exterior  grid  spacing. 


T able  2.  Elastic  properties  used  in  the  Shoal  calcu 


Density 

2600  kg/m^ 

P-veloeity 

5175  tn/s 

S-veloeity 

3026  m/s 

ation 


Figure  14  shows  time  domain  RVP  and  souree  speetra  calculated  in  two  ways,  first  with  a  ID 
spherically  symmetric  code  and  second  with  CRAM3D  derived  from  the  first  second  of  a  radial 
waveform  at  shot  level.  The  results  are  similar,  but  the  3D  eode  breaks  the  perfect  symmetry  of 
the  spherically  symmetric  calculation,  causing  reduced  rebound  and  extending  motion  over  a 
longer  time  period.  This  eauses  a  less  peaked  spectrum  and  higher  statie  level. 


Figure  14.  Source  function  for  material  properties  used  in  the  Shoal  ealeulation  ealeulated  using 
a  ID  spherieally  symmetric  code  and  using  a  shot  level  free-field  waveform  from  Cram3D.  Left: 
speetral  amplitude;  Right:  RVP  time  series. 


4.2.1  Shoal:  Near  Field  Waveforms 

Figure  15  shows  the  observed  shot  level  radial  ground  motion  at  the  three  stations  deseribed 
above  and  shown  in  Figure  12.  The  motion  is  very  different  at  the  three  stations.  In  partieular 
station  PM-3  had  a  large  positive  final  displacement,  while  at  the  other  two  stations  the  final 
displaeement  was  small,  and  in  faet  slightly  negative  at  PM-1. 
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Shoal  -  Radial  Velocity 


Shoal  -  Radial  Displacement 


Figure  15.  Observed  near  field  (shot  level)  velocity  (left)  and  displacement  (right)  at  the  three 
Shoal  stations. 

Figure  16  shows  the  displacement  waveforms  observed  at  the  three  stations  compared  with 
calculated  waveforms  at  the  same  depth  in  the  directions  of  minimum  and  maximum  horizontal 
stress.  As  noted  above,  one  of  the  surprising  aspects  of  the  data  from  this  event  is  the  strong 
directionality  of  the  shot-level  seismic  waveforms,  in  particular  the  fact  that  the  final 
displacements  were  so  different  at  the  three  stations.  We  see  a  similar  effect  in  our  calculations. 
The  late  time  displacement  is  amplified  in  the  direction  of  minimum  principal  stress  compared  to 
the  other  directions  (smaller  than  the  observations,  but  otherwise  very  consistent).  PM-3,  the 
station  with  the  largest  displacement,  is  in  the  direction  of  the  regional  tectonic  minimum 
principal  stress,  so  prestress  is  very  likely  a  strong  contributor  to  the  differences  in  the  observed 
displacements. 


Approved  for  public  release;  distribution  is  unlimited. 
16 


Shoal  -  Radial  Displacement 
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Figure  16.  Left:  observed  displaeement  waveforms  at  the  three  shot-level  Shoal  reeording 
stations.  Right:  ealeulations  with  prestress  in  the  direetions  of  maximum  (Az=0)  and  minimum 
(Az=90)  horizontal  stress.  Both  data  and  ealeulations  show  enhaneed  late  time  displaeements  in 
the  direetion  of  minimum  prineipal  stress. 

4.2.2  Shoal:  Regional  Waveforms  and  Surface  Waves 

We  calculated  regional  waveforms  using  full  waveform  Green’s  functions  at  250  km  distance 
from  two  calculations,  the  first  with  no  tectonic  stress  and  the  second  with  tectonic  stress  as 
described  above.  In  the  calculation  without  tectonic  release,  the  seismic  radiation  is  isotropic  and 
there  is  no  motion  on  the  tangential  component  (Figure  17).  However,  with  tectonic  release  a 
strong  Love  wave  is  generated  and  the  Rayleigh  waves  have  a  radiation  pattern,  with  amplitudes 
increased  in  the  direction  of  extensional  stress  and  decreased  in  the  direction  of  (smaller) 
compressive  stress  (Figure  18).  Note  that  although  tectonic  release  causes  variability  in  the 
Rayleigh  waves  and  tangential  motion  not  present  at  all  without  it,  the  regional  body  waves  are 
not  significantly  altered  by  tectonic  release. 
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Figure  17.  CRAM3D  calculation  of  Shoal  without  prestress  propagated  to  250  km  using  the 
representation  theorem.  These  are  broadband  full  waveform  seismograms  calculated  using 
wavenumber  integration.  Seismic  waveforms  are  independent  of  observation  direction. 
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Figure  18.  CRAM3D  calculation  of  Shoal  with  tectonic  prestress  propagated  to  250  km  using  the 
representation  theorem.  Now  there  is  a  love  wave  and  variability  in  the  Rayleigh  wave  with 
azimuth.  Zero  degrees  is  the  direction  of  extension. 
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4.2.3  Shoal:  Long  Period  Surface  Waves 


To  see  the  long  period  behavior  more  clearly,  we  generated  seismograms  using  the 
representation  theorem  with  fundamental  mode  surface  wave  Green’s  functions.  A  problem 
occurs  when  trying  to  apply  the  representation  theorem  to  long  period  surface  waves.  The 
Green’s  function  for  a  vertical  force  becomes  large  at  long  periods,  so  any  small  vertical  motion 
left  in  the  grid  is  amplified  and  dominates  the  signal.  In  fact,  this  term  must  vanish  in  the  long 
period  limit.  To  resolve  this  problem  we  extend  the  time  series  by  0.1  seconds  and  adjust  the 
stresses  on  the  side  and  bottom  so  that  the  total  force  and  impulse  on  the  monitoring  surface  is 
zero  at  the  end  of  the  calculation.  This  problem  was  described  for  axisymmetric  problems  and  a 
similar  technique  applied  by  Bache  et  al.  (1982). 

Fundamental  mode  surface  waves  calculated  at  250  km  and  low  pass  filtered  at  0.1  Hz  are  shown 
in  Figure  19.  The  Rayleigh  wave  is  increased  in  amplitude  by  about  60%  in  the  direction  of 
extension  compared  to  the  direction  of  (mild)  compression.  A  Love  wave  is  generated  with 
maximum  amplitude  about  half  of  the  Rayleigh  wave  amplitude.  We  compared  the  Rayleigh  and 
Love  waves  with  observations  of  surface  waves  by  Lambert,  Flynn  and  Archambeau  (1972)  and 
found  that  the  calculations  are  consistent  with  the  observed  radiation  pattern  (Figure  20). 
Although  the  observed  effect  appears  to  be  larger  than  the  calculated  effect,  the  radiation  pattern 
clearly  shows  amplification  in  the  direction  of  extension  and  a  Love  wave  radiation  pattern 
consistent  with  the  Rayleigh  wave  radiation  pattern. 


Vertical  Velocity  Tangential  Velocity 


Figure  19.  CRAM3D  calculation  of  Shoal  with  prestress  propagated  to  250  km  using 
fundamental  mode  Rayleigh  and  Love  wave  Green’s  functions  low  pass  filtered  at  0.1  Hz. 
Rayleigh  waves  vary  in  amplitude  with  angle  and  a  Love  wave  is  generated. 
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Figure  20.  From  Lambert,  Flynn  and  Archambeau  (1972).  The  Rayleigh  and  Love  wave 
radiation  pattern  observed  from  the  Shoal  explosion  is  consistent  with  the  calculation. 


We  can  estimate  the  size  of  the  calculated  tectonic  source  by  comparing  the  surface  waves  for 
the  calculations  with  and  without  prestress.  We  find  that  a  moment  tensor  with  Mxx=7.5xl0'"^ 
and  Myy=-2.5xl0*'^  added  to  the  calculation  without  prestress  matches  the  Love  wave  and  the 
Rayleigh  wave  variation  for  the  calculation  with  prestress  very  well  (Figure  21).  The  isotropic 
moment  for  the  ID  calculation  in  the  same  material  (Figure  14)  had  a  moment  of  1.9x10*^,  so  the 
perturbation  in  the  horizontal  components  is  about  half  of  the  isotropic  moment. 
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Figure  21.  Left:  Fundamental  mode  Rayleigh  and  Love  wave  spectra  from  the  Shoal  calculation 
with  prestress  at  250  km.  Right:  Calculation  without  prestress  plus  point  moment  tensor  solution 
as  discussed  in  the  text. 
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We  can  take  this  analysis  a  bit  further  and  determine  the  effect  of  the  free  surface  on  the  vertical 
and  horizontal  moment  tensor  components.  Figure  22  shows  a  comparison  of  surface  wave 
spectra  from  the  3D  calculation  compared  with  spectra  from  several  ID  calculations.  The 
isotropic  moment  of  the  3D  source  is  well  defined  by  the  long  period  limit  of  the  surface  wave 
spectrum  and  is  equal  to  2.6x10*^  N-m.  As  discussed  earlier,  there  is  more  overshoot  in  the  ID 
than  the  3D  source  which  reduces  the  low  frequency  level  (Figure  14),  so  we  multiply  the 
spectrum  from  the  point  explosion  source  convolved  with  the  ID  source  function  by  1.4  to  make 
the  long  period  levels  the  same.  At  frequencies  below  0.1  Hz,  all  of  the  spectra  are 
indistinguishable  and  consistent  with  an  isotropic  point  source  solution.  At  higher  frequencies 
(-0.4  Hz)  there  is  a  divergence  of  the  spectra  for  two  reasons:  in  the  ID  case  it  is  because  of  the 
spectral  shape  caused  by  rebound;  in  the  3D  case  it  is  because  of  the  free  surface  interaction.  The 
spectra  can  be  matched  by  increasing  the  Mzz  component  and  also  increasing  the  Mxx  and  Myy 
components  by  a  lesser  amount.  The  red  curve  in  the  figure,  which  fits  the  full  spectral  shape 
quite  well,  has  Mzz  increased  by  a  factor  of  2.5  over  the  apparent  isotropic  moment,  while  Mxx 
and  Myy  are  increased  by  a  factor  of  1.5,  so  the  final  moment  tensor  is  Mzz=6.4xl0*^, 
Mxx=Myy=3.8xl0'^.  The  isotropic  moment  increases  to  4.67x10^^  and  there  is  a  CLVD 
component  caused  by  the  difference  between  the  horizontal  and  vertical  stresses.  This  is  similar 
to  the  “damage”  mechanism  that  has  been  discussed  by  Patton  and  Taylor  (2011).  Note, 
however,  that  although  the  CLVD  caused  by  the  free  surface  interaction  reduces  the  surface 
wave  amplitude,  the  isotropic  moment  is  also  increased  and  so  the  net  effect  on  the  surface  wave 
amplitude  is  small.  The  surface  wave  generated  by  the  ID  source  is  actually  smaller  by  a  factor 
of  1 .4  than  the  3D  solution  as  discussed  above,  but  a  direct  comparison  is  not  possible  because 
the  overshoot  in  the  ID  case  causes  the  surface  wave  amplitude  to  be  lower,  and  that  effect  is 
reduced  in  the  3D  case. 


Shoal  Surface  Waves  3D+Mxx+Myy+Fac*Mzz 


Figure  22.  Comparison  of  ID  and  3D  surface  wave  spectra.  All  are  fundamental  mode  spectra  at 
a  distance  of  250  km.  The  black  dashed  line  labeled  “3D”  is  the  spectrum  from  the  3D 
calculation  without  tectonic  release.  The  black  dotted  line  is  from  a  ID  calculation  convolved 
with  the  spherically  symmetric  source  function  multiplied  by  1.4.  The  three  colored  curves  are 
ID  point  source  calculations  constrained  to  have  the  same  long  period  limit  as  the  3D 
calculation,  but  with  different  values  of  Mzz  vs.  Mxx  and  Myy  as  discussed  in  the  text. 
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4.2.4  Shoal:  Far  Field  Body  Waves 


We  calculated  far-field  body  waves  from  the  Shoal  calculations  with  and  without  tectonic  release 
to  examine  the  effect  on  body  wave  signatures  and  generation  of  S  waves.  The  results  show  only 
small  differences  in  the  P  and  SV  waves  (SV  waves  are  primarily  generated  by  P  to  S  conversion 
at  the  free  surface),  however  tectonic  release  generates  a  far-field  SH  wave  that  is  absent  in  the 
case  without  tectonic  release.  Figure  23,  Figure  24  and  Figure  25  show  P,  SV  and  SH  waves 
calculated  in  the  North,  Northeast  and  East  directions,  where  North  is  the  direction  of  extension 
in  the  calculation.  All  waveforms  were  calculated  at  a  distance  of  1000  km  in  a  granite  half  space 
and  were  filtered  with  a  causal  low  pass  filter  at  5  Hz.  All  waveforms  are  downgoing  waves  as  if 
they  were  measured  1000  km  below  the  source  in  a  perfect  granite  medium  undisturbed  by 
reflections  (other  than  above  the  source)  or  mantle  complexity,  so  these  waveforms  provide  a 
very  clear  image  of  the  seismic  source. 
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Figure  23.  P  waves  from  Shoal  calculation  with  tectonic  release  (black)  and  without  (red).  There 
are  very  small  differences  in  the  waveforms.  Units  are  velocity  normalized  to  4x10'"^  m/s.  TOA  is 
takeoff  angle  from  the  vertical.  N,  NE  and  E  indicate  the  direction  to  the  receiver  as  discussed  in 
the  text. 


Approved  for  public  release;  distribution  is  unlimited. 
22 


SV:  10  TOA,  N 


SV:  20  TOA,  N 


SV:  30  TOA,  N 


SV:  10TOA,  NE 


-1 

1.5  2  2.5  3  3.5 

SV;  10  TOA,  E 


-1 


1.5  2  2.5  3  3.5 

SV:  20  TOA,  NE 


SV:  20  TOA,  E 


SV:  30  TOA,  NE 


SV:  30  TOA,  E 


1.5  2  2.5  3  3.5 


Figure  24.  SV  waves  from  Shoal  calculation  with  tectonic  release  (black)  and  without  (red). 
There  are  very  small  differences  in  the  waveforms.  Units  are  velocity  normalized  to  4x10"^  m/s. 
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Figure  25.  SH  waves  from  Shoal  calculation  with  tectonic  release  (black)  and  without  (red).  In 
the  NE  direction,  a  substantial  SH  wave  is  generated  in  the  calculation  with  tectonic  release  that 
is  not  present  in  the  calculation  without  tectonic  release.  Again,  units  are  velocity  normalized  to 
4x10“^  m/s,  but  note  that  the  limits  are  different  from  the  previous  figures.  The  SH  wave  in  the 
NE  direction  is  about  10%  the  size  of  the  SV  waves.  SH  is  zero  in  the  north  and  east  directions. 

Although  the  effect  of  tectonic  release  on  body  wave  amplitudes  appears  to  be  small,  the  tectonic 
release  signal  in  the  body  wave  can  be  identified  by  subtracting  the  body  waves  from  the  two 
calculations.  Figure  26  shows  the  far  field  P  waves  along  the  directions  of  extension  and 
compression.  As  expected,  the  extension  causes  an  increase  in  signal.  The  tectonic  signal  is 
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about  10%  the  size  of  the  direet  signal  and  delayed  slightly  in  time,  so  the  peak  amplitude  of  the 
signal  is  inereased  only  very  slightly.  In  the  direetion  of  compression  (recall  that  compression 
was  considerably  weaker  than  extension  in  the  other  direction)  the  effect  is  smaller  with  the  main 
effect  being  a  slight  amplitude  decrease  at  about  1 .9  seconds. 
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Figure  26.  Top  row:  Far  field  P  waves  in  the  X  and  Y  directions  -  X  (North)  is  the  direction  of 
extension.  Bottom  row:  Velocity  differences  subtracting  the  waveform  without  prestress  from  the 
waveform  with  prestress.  Both  are  for  a  30  degree  takeoff  angle.  Units  are  velocity  normalized  to 
4x10"^  m/s. 


The  effect  is  similar  but  larger  for  far  field  SV  (Figure  27),  which  shows  that  the  SV  waveforms 
from  tectonic  release  are  about  20%  of  the  direct  SV.  It  also  starts  earlier,  because  the  “direct” 
SV  wave  is  dominated  by  pS  and  so  is  delayed  by  the  travel  time  to  the  surface,  while  the 
tectonic  SV  is  coincident  with  the  explosion.  Figure  28  shows  that  an  SH  wave  is  generated  in 
the  calculation  with  tectonic  release  that  is  not  present  in  the  calculation  without  tectonic  release. 
The  SH  wave  at  its  maximum  at  45  degrees  between  the  Shmin  and  SHmax  direction  is  about 
10%  the  size  of  the  SV  waves.  SH  is  zero  in  the  direction  of  principal  stresses  which  are  axes  of 
symmetry. 
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Figure  27.  Top  row:  Far  field  SV  waves  in  the  X  and  Y  directions  -  X  (North)  is  the  direction  of 
extension.  Bottom  row:  Velocity  differences  subtracting  the  waveform  without  prestress  from  the 
waveform  with  prestress.  Both  are  for  a  30  degree  takeoff  angle.  Units  are  velocity  normalized  to 
4x10"^  m/s. 
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Figure  28.  Far  field  SH-waves  calculated  at  10  (left),  20  (middle)  and  30  (right)  degree  take-off 
angles  from  Shoal  calculations  with  and  without  prestress;  45  degrees  from  the  direction  of 
extension. 


In  Figure  29,  we  compare  the  far  field  P-waves  calculated  for  the  explosion  without  tectonic 
release  with  a  calculation  for  a  spherically  symmetric  explosion  source  calculated  in  an  infinite 
medium  with  the  same  material  properties  as  the  3D  calculation,  and  placed  at  the  same  depth  in 
the  same  structure.  The  difference  between  these  waveforms  corresponds  to  the  nonlinear  effects 
caused  by  the  free  surface  and  the  variation  of  overburden  pressure  with  depth.  These  effects 
include  spall,  material  damage  and  nonlinear  deformation.  The  main  difference  is  a  reduction  in 
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the  pP  amplitude  and  a  longer  duration  in  the  3D  case.  Similar  results  were  found  for  two- 
dimensional  calculations  by  Stevens  et  al  (1991). 
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Figure  29.  Far  field  P-waves  calculated  for  a  spherically  symmetric  explosion  source  (red) 
compared  with  far  field  P-waves  from  the  Shoal  calculation  without  tectonic  release  (black). 
Both  calculations  used  the  same  material  properties  and  earth  structure.  From  left  to  right,  take¬ 
off  angles  of  10,  20  and  30  degrees. 


Figure  30  shows  the  same  comparison  for  SV  waves.  The  3D  calculation  shows  larger  SV 
waves,  particularly  at  shallower  takeoff  angles,  and  a  precursor  not  present  in  the  spherically 
symmetric  case.  The  largest  part  of  the  SV  wave  is  again  the  pS  phase,  modified  by  the  nonlinear 
interaction  with  the  free  surface.  The  precursor  is  a  direct  S  wave  generated  by  the  explosion  in 
the  3D  medium. 
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Figure  30.  Far  field  SV-waves  calculated  for  a  spherically  symmetric  explosion  source  (red) 
compared  with  far  field  SV-waves  from  the  Shoal  calculation  without  tectonic  release  (black). 
Both  calculations  used  the  same  material  properties  and  earth  structure.  From  left  to  right,  take¬ 
off  angles  of  10,  20  and  30  degrees. 
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4.2.5  Shoal:  Hypothetical  Compressive  Prestress 


Shoal  is  in  a  region  of  extension  whieh,  as  shown  above,  increases  surface  wave  amplitudes, 
although  the  effect  is  not  large  with  the  smaller  compressive  prestress  in  the  other  direction.  In 
this  section  we  consider  the  effects  that  would  occur  if  the  same  explosion  in  the  same  medium 
were  conducted  in  a  region  in  compressive  prestress.  In  this  calculation,  we  left  the  secondary 
prestress  the  same  and  changed  the  dominant  prestress  from  tensile  to  compressive  (also  rotated 
90”),  so  the  medium  is  balanced  against  friction  for  a  thrust  fault  instead  of  a  normal  fault.  The 
initial  stress  state  for  the  calculation  is  shown  in  Figure  31. 


Figure  31.  Stress  state  for  calculation  of  Shoal  with  compressive  prestress.  The  dashed  line  is  the 
critical  SHmax  stress  limited  by  friction. 

The  broadband  regional  waveforms  calculated  from  this  simulation  are  shown  in  Figure  32.  At 
first  glance  they  look  similar  to  the  case  with  tensile  prestress  -  there  is  a  clear  Love  wave  and 
the  Rayleigh  wave  varies  with  azimuth  while  the  body  waves  do  not,  and  in  fact  the  entire 
wavetrain  ahead  of  the  surface  waves  is  very  similar  to  the  other  cases.  Figure  33  shows  the 
fundamental  mode  waveforms  lowpass  filtered  at  0. 1  Hz  and  compared  with  the  same  waveform 
for  the  tensile  prestress  case.  The  variation  in  amplitude  is  opposite  of  the  tensile  case.  As 
expected,  it  is  smallest  in  the  direction  of  maximum  compressive  stress  (90'’). 

Figure  34  -  Figure  36  show  far-field  body  waves  from  the  compressive  calculation  compared 
with  the  other  calculations.  The  results  are  very  similar  to  the  tensile  prestress  calculation.  There 
is  little  effect  on  the  initial  waveform,  but  the  later  part  of  the  waveform  that  is  affected  by  the 
free  surface  interaction  is  delayed  and  modified  by  prestress.  The  effects  described  earlier  are 
somewhat  larger  for  compressive  prestress  compared  to  tensile  prestress. 
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Figure  32.  Broadband  waveforms  from  the  Shoal  ealculation  with  compressive  prestress. 


Vertical  Velocity  Tangential  Velocity 


/V 

-1  0Deg| 

-2 

xIO^  - - - - — 

A 

ynxr^ 

.  1  45  Deg 

Ny  W'  ^ 

-2 

.... 

.... 

X  10 

2 

A_  -  0 

X  10 

.  1  90  Deg 

-2 

xIO^  ^ — 

X  10  ^^  ‘ - ^ ^ ^ — 

Figure  33.  Fundamental  mode  waveforms  at  250  km  from  the  calculation  with  compressive 
prestress  low  pass  filtered  at  0.1  Hz.  The  bottom  row  shows  the  waveforms  from  the  Shoal 
calculation  described  in  the  previous  section  degrees  for  comparison.  Surface  wave  amplitude  for 
the  compressive  case  in  the  direction  of  compressive  prestress  is  reduced. 
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Figure  34.  Top  row:  Far  field  P  waves  in  the  X  and  Y  directions  -  X  (North)  is  the  direction  of 
extension.  Bottom  row:  Velocity  differences  subtracting  the  waveform  without  prestress  from  the 
waveforms  with  the  Shoal  prestress  and  hypothetical  compressive  prestress.  Both  are  for  a  30 
degree  takeoff  angle.  Units  are  velocity  normalized  to  4x10"^  m/s. 
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Figure  35.  Top  row:  Far  field  SV  waves  in  the  X  and  Y  directions  -  X  (North)  is  the  direction  of 
extension.  Bottom  row:  Velocity  differences  subtracting  the  waveforms  without  prestress  from 
the  waveforms  with  the  Shoal  prestress  and  hypothetical  compressive  prestress.  Both  are  for  a  30 
degree  takeoff  angle.  Units  are  velocity  normalized  to  4x10"^  m/s. 
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Figure  36.  Far  field  SH-waves  calculated  at  10  (left),  20  (middle)  and  30  (right)  degree  take-off 
angles  from  Shoal  calculations  without  prestress  and  the  waveforms  with  the  Shoal  prestress  and 
hypothetical  compressive  prestress;  45  degrees  from  the  direction  of  extension. 


4.3  3D  Earth  Heterogeneity 

All  explosions  generate  shear  waves,  including  even  overburied  and  decoupled  explosions.  A 
possible  explanation  for  this  is  that  the  extra  degree  of  freedom  provided  in  three  dimensions 
makes  generation  of  shear  waves  much  easier  than  with  the  restriction  to  axisymmetric 
geometry.  In  other  words,  spherical  and  axisymmetric  calculations  artificially  suppress  the  shear 
waves.  In  this  section  we  try  to  estimate  the  magnitude  of  this  effect  by  incorporating 
heterogeneous  material  properties  into  3D  nonlinear  calculations. 

The  spatial  distribution  of  the  material  properties  of  the  earth’s  crust  is  composed  of  step-like 
changes  and  short  wavelength  random  fluctuations,  which  strongly  affect  the  propagation  of 
elastic  waves.  The  most  difficult  part  is  defining  the  characteristics  of  the  heterogeneous  media. 
There  is  very  little  information  available  on  strength  heterogeneity,  for  example.  There  is 
considerably  more  information  available  on  velocity  and  density  heterogeneity,  primarily  from 
well  logs  and  core  samples.  We  are  using  these  in  two  ways:  1)  to  develop  velocity  and  density 
models  based  on  the  observed  distributions;  and  2)  to  develop  a  strength  distibution  model  under 
the  assumption  that  the  strength  heterogeneity  is  similar  to  the  velocity  distribution. 

Nine  3D  finite  element  calculations  were  performed  of  the  Shoal  explosion  with  CRAM3D  to 
investigate  body  and  surface  wave  generation  with  heterogeneous  material  properties  in  the  shot 
region  (Table  3).  Three  pairs  of  the  finite  element  calculations  used  distributions  of  density,  bulk 
modulus,  and  shear  modulus  with  statistical  characteristics  the  same  as  those  derived  from  well 
YT2  log  data  from  Kyushu,  Japan  (Figure  37).  These  calculations  are  referred  to  as  the 
“lithosphere”  and  “weakened  lithosphere”  calculations  below  since  the  material  distributions  are 
modeled  after  actual  measurements.  To  test  the  sensitivity  of  the  waveforms  to  the  spatial 
distribution  function,  two  calculations  were  performed  with  slightly  modified  spatial 
distributions  (“modified  lithosphere”). 
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Figure  37.  Taken  from  Shiomi  et  al.  1997:  Well  logs  showing  P-  and  S-wave  velocities  and  mass 
density  versus  depth  for  well  YT2  in  Kyushu,  Japan  for  depths  of  600  m  to  1710  m. 


Table  3.  Nine  3D  finite  element  calculations  were  performed  to  investigate  3D  Heterogeneity. 


Calculation 

Heterogeneities 

Strength 

weakening 

Number  of 
Calcs. 

Uniform  media 

None 

No 

1 

Lithosphere 

As  in  well  YT2 

No 

3 

Lithosphere  (weakened) 

As  in  well  YT2 

Yes 

3 

Modified  lithosphere 

Modified  from  YT2 

No 

1 

Mod.  lith.  (weakened) 

Modified  from  YT2 

Yes 

1 

Bulk  modulus  (GPa) 
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Figure  38.  Top  panels:  Examples  of  the  randomly  generated  material  property  distributions  for 
one  set  of  the  “lithosphere”  calculations  for  a  2  km  x  2  km  horizontal  slice  at  100  m  depth.  The 
resolution  (cell  size)  of  the  material  property  distribution  is  20  m.  The  statistical  properties  of 
these  distributions  are  identical  to  those  in  well  YT2.  Bottom  panels.  Same  as  the  top  panels,  but 
for  the  “modified  lithosphere”  calculations. 
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Although  the  YT2  well  log  data  only  reveal  the  distribution  of  material  properties  with  depth,  we 
assume  that  the  same  autocorrelation  function  applies  in  all  spatial  dimensions.  Figure  38  shows 
horizontal  slices  of  the  material  properties  from  two  of  the  calculations.  The  top  and  bottom 
panels  correspond  to  the  spatial  distribution  from  the  YT2  well  log  data,  and  from  the  slightly 
modified  distribution,  respectively.  The  bulk  modulus  spans  10-58  GPa  with  a  median  value  of 
27.4  GPa,  the  shear  modulus  spans  6-25  GPa  with  a  median  value  of  14.4  GPa,  and  density  spans 
1800-2600  kg/m^  with  a  median  value  of  2370  kg/m^  (Figure  39).  The  distributions  of  bulk  and 
shear  modulus  are  highly  correlated  with  one  another  (0.92  correlation  coefficient),  and  the 
density  distribution  is  correlated  to  the  moduli  but  to  a  lesser  extent  (0.3 1  correlation  coefficient). 
Because  of  this  correlation,  many  of  local  maxima  in  the  bulk  modulus  distribution,  for  example, 
align  with  the  maxima  of  the  shear  modulus  distribution. 


Shear  Modulus  (GPa)  Bulk  Modulus  (GPa)  Density  (kg/m^) 


Figure  39.  The  distributions  of  the  properties  derived  from  the  YT2  well  log  data. 

The  autocorrelation  functions  of  the  well  log  data  and  the  two  simulated  distributions  are  shown 
in  Figure  40.  The  bulk  modulus  and  shear  modulus  exhibit  very  similar  autocorrelation 
functions;  they  are  spatially  correlated  over  distances  of  ~80  m.  The  density  fluctuations  are 
correlated  over  distances  of  ~30  m.  For  the  modified  spatial  distribution  functions,  the  bulk  and 
shear  moduli  are  correlated  over  slightly  larger  distances,  and  the  density  is  correlated  over 
slightly  smaller  distances  relative  to  the  other  calculations. 


Figure  40.  Left  panel:  Autocorrelation  function  of  P  wave  velocity,  S  wave  velocity,  and  density 
of  the  data  from  well  YT2.  The  numerals  are  the  mean-squared  fractional  fluctuation  for  each 
property.  Middle  panel:  Ensemble  average  of  the  autocorrelation  functions  of  bulk  modulus, 
shear  modulus,  and  density  used  in  the  “lithosphere”  calculations.  Right  panel:  Same  as  the 
middle  panel,  but  for  the  “modified  lithosphere”  calculation. 
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In  all  the  calculations,  the  heterogeneities  are  tapered  to  the  mean  values  between  600  m  and  800 
m  from  the  shot  point  along  x  and  y  dimensions,  and  between  200  m  and  400  m  below  the  shot 
point  in  the  z  dimension  (heterogeneities  extend  to  the  surface).  The  purpose  of  the  tapering  is  to 
have  constant  material  properties  along  the  monitoring  surface.  The  monitoring  surface  is  where 
the  displacements  and  stresses  are  saved  to  propagate  the  wavefield  from  the  boundary  of  the 
calculation  to  an  arbitrary  exterior  observation  point  using  the  elastodynamic  representation 
theorem. 

In  addition  to  investigating  the  body  waves  and  surface  waves  generated  from  two  spatial 
correlation  functions,  we  also  tested  the  effects  of  these  distributions  with  the  inclusion  of 
correspondent  weakening  of  the  material  properties.  Material  weakening  (shock  damage)  is 
incorporated  through  a  reduction  in  both  strength  and  shear  modulus  based  on  the  maximum 
shear  strain  seen  by  the  rock  element.  Although  material  weakening  is  always  included  in  the 
calculations  once  a  minimum  amount  of  plastic  work  has  been  surpassed,  for  the  “weakened” 
calculations  the  reduced  shear  modulus  and  the  shear  strain  thresholds  that  control  the  weakening 
are  lowered.  Specifically,  for  each  cell  five  parameters  *  are  all  scaled  by  the  ratio  of  the  shear 
modulus  in  the  cell  to  the  maximum  shear  modulus  used  in  the  calculation  (25  GPa).  (See 
Section  1.5.3  of  the  CRAM3D  Users  Manual  for  a  more  detailed  explanation  of  the  shock 
damage  parameters.) 

Finally,  we  also  ran  a  calculation  with  uniform  material  properties  in  order  to  remove  the  effects 
of  the  heterogeneities  on  the  waveforms.  The  mean  values  of  the  varied  parameters  were  used  in 
this  calculation. 

Table  3  provides  a  summary  of  the  nine  3D  finite  element  calculations.  For  each  calculation, 
body  waves  were  calculated  for  take-off  angles  of  10,  20,  and  30  degrees  along  4  azimuths  (to 
the  north,  east,  south,  and  west).  Surface  waves  were  calculated  for  4  azimuths. 

4.3.1  Heterogeneous  Source  Region:  Far  Field  Body  Waves 

The  body  waveforms  for  the  “lithosphere”  and  the  “weakened  lithosphere”  calculations  are 
shown  in  Figure  41.  In  addition,  we  show  the  “lithosphere”  body  waves  with  the  “uniform”  body 
waves  subtracted  to  try  to  isolate  the  effects  of  the  heterogeneities.  Root-mean-squared  (RMS) 
amplitudes  of  the  body  waves  were  calculated  using  2  seconds  of  data  containing  the  entire 
waveforms  (Table  4). 


^  deinit,  defull,  epinit,  epfull,  gdam;  See  the  CRAM3D  Users  Manual. 
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Figure  41.  Top:  Body  waves  for  a  10  degree  take-off  angle,  low-pass  filtered  at  20  Hz,  for  4 
directions  for  each  of  the  3  “lithosphere”  calculations.  The  amplitudes  of  the  waveforms  are 
scaled  by  le-4  and  are  for  1000  km  range.  Time  is  on  the  x-axis;  a  negative  offset  of  0.333  s  has 
been  applied  to  the  shear  waves  to  align  the  start  of  the  motion.  Middle:  Same  as  the  top  panel, 
but  for  the  3  “weakened  lithosphere”  calculations.  Bottom:  Same  as  the  top  panels,  but  with  the 
“uniform”  media  calculation  subtracted  from  the  waveforms. 

Including  heterogeneity  in  the  material  properties  of  the  source  region  has  only  a  subtle  effect  on 
the  RMS  amplitudes  of  P  and  SV  waves  (1-6  percent  increases  relative  to  uniform  media).  SH 
waves,  which  are  non-existent  in  uniform  media,  are  generated  with  RMS  amplitudes  about  20- 
33  percent  as  large  as  SV  waves.  When  concurrent  weakening  of  the  material  in  also  included, 
the  RMS  amplitudes  for  all  body  wave  types  are  increased  on  average.  P  wave  amplitudes 
increase  10-20  percent  relative  to  uniform  and  heterogeneous  cases  without  weakening.  SH 
amplitudes  increase  substantially,  by  70-150  percent.  The  RMS  amplitudes  of  SV  waves  also 
increase  nearly  70-80  percent  for  a  10  degree  take-off  angle,  but  the  amplification  from  material 
weakening  is  greatly  diminished  at  larger  take  off  angles,  with  only  about  10  percent  increase  at 
a  30  degree  take-off  angle. 
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Table  4.  R] 

MS  amp 

itudes  0] 

f  body  waves  for  the  nine  CRAM3D  calculations. 

Calculation 

PIO 

P20 

P30 

SHIO 

SH20 

SH30 

SVIO 

SV20 

SV30 

Uniform 

4.19 

4.09 

3.98 

0 

0 

0 

6.53 

8.46 

7.58 

Lithosphere  1 

4.41 

4.19 

3.98 

1.37 

1.16 

1.21 

6.70 

8.41 

7.54 

Lithosphere  2 

4.43 

4.16 

3.97 

1.36 

1.17 

1.22 

6.76 

8.42 

7.55 

Lithosphere  3 

4.26 

4.07 

3.88 

1.28 

1.44 

1.23 

6.88 

8.24 

7.21 

Modified  lith.  1 

4.32 

4.17 

4.00 

1.77 

1.62 

1.59 

5.96 

7.93 

7.32 

Lith.  1  (weakened) 

5.13 

4.74 

4.42 

2.34 

2.82 

2.11 

12.23 

9.67 

8.45 

Lith.  2  (weakened) 

5.09 

4.70 

4.42 

2.33 

2.81 

2.09 

12.19 

9.68 

8.46 

Lith.  3  (weakened) 

4.76 

4.82 

4.62 

3.19 

2.48 

2.24 

11.87 

9.90 

8.49 

Mod.  lith.  1  (wknd) 

4.95 

4.70 

4.48 

3.10 

3.62 

2.70 

10.70 

9.27 

8.15 

Note. — RMS  amplitudes  for  P,  SH,  and  SV  waves  at  10,  20,  and  30  degree  take-off  angles.  Eaeh  measurement  is 
derived  from  the  average  of  4  azimuths  (north,  east,  south,  west).  Amplitudes  are  sealed  by  le-5. 


Because  P  and  SV  waves  are  present  even  in  the  absence  of  heterogeneity,  only  subtle 
differences  with  azimuth  are  observed.  The  largest  differences  in  the  waveforms  occur  where  the 
gradient  of  the  waveform  is  largest  (Figure  41,  bottom),  indicating  that  the  waves  are  primarily 
affected  by  propagation  speed  fluctuations.  Similarly,  the  particular  random  instance  drawn  from 
the  “lithosphere”  or  “modified  lithosphere”  distribution  functions  (1,  2,  or  3)  has  a  relatively 
minor  effect  on  the  P  and  SV  waves.  On  the  contrary,  SH  waves  show  ~20-30%  variation  in 
RMS  amplitude  with  azimuth. 


Table  5.  Azimuthal  variation  (percent)  of  body  waves  in  the  presence  of  heterogeneity  (la) 


Calculation 

PIO 

P20 

P30 

SHlO 

SH20 

SH30 

svio 

SV20 

SV30 

Lithosphere 

2 

2 

3 

13 

19 

21 

6 

4 

4 

Modified  lith. 

1 

3 

5 

18 

16 

27 

7 

8 

4 

Lith.  (weakened) 

3 

2 

3 

25 

17 

18 

7 

9 

3 

Mod.  lith.  (wknd) 

1 

3 

3 

33 

11 

5 

3 

3 

2 

Note. — Tabulation  of  the  percent  variation  with  azimuth  combines  the  3  “lithosphere”  and  “weakened  lithosphere” 
calculations  (12  azimuths). 


The  power  spectral  density  of  the  waveforms  for  P  and  SV  waves  are  only  slightly  different  in 
the  “weakened”  calculations  relative  to  the  non-weakening  calculations  (Figure  42).  Both  types 
of  waveforms  show  a  subtle  shift  of  energy  from  high  frequencies  to  lower  frequencies.  SH 
waves,  on  the  other  hand,  show  a  more  distinct  shift  in  energy  to  longer  periods. 
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P:  10TOA 


P:  20  TO  A 


P:  30  TO  A 


Figure  42.  Normalized  power  spectral  density  of  the  body  waves  for  the  4  types  of  scenarios 
under  investigation.  Distinct  wave  types  (take-off  angles)  are  shown  in  each  row  (column). 


4.3.2  Heterogeneous  Source  Region:  Broadband  Waveforms  and  Surface 
Waves 


Broadband  regional  waveforms  at  250  km  were  calculated  from  the  nine  CRAM3D  calculations. 
As  seen  in  Table  6,  heterogeneous  material  properties  in  the  source  region  (without  weakening) 
generate  Love  waves,  but  have  only  a  small  effect  on  the  Rayleigh  waves.  However,  the  RMS 
amplitudes  of  the  Love  waves  are  ~3x  greater  when  enhanced  weakening  of  the  material 
properties  is  included  in  the  calculations.  Rayleigh  wave  amplitudes  increase  a  more  modest 
~40%  with  weakening.  Comparison  of  Figure  43  and  Figure  44  clearly  reveals  the  significant 
increase  in  surface  wave  amplitudes  with  material  weakening.  Similar  to  SH  waves,  the  standard 
deviation  of  the  RMS  amplitudes  of  the  Love  waves  over  all  azimuths  is  about  20%,  independent 
of  whether  or  not  enhanced  material  weakening  in  included  in  the  calculation.  The  variation  with 
azimuth  of  the  Rayleigh  waves  is  very  small  (~3%)  in  all  heterogeneity  cases. 
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Table  6.  Surface  waves  with  uniform  and  heterogeneous  material  properties  in  the  source  region. 


Calculation 

Love  waves  RMS 
amplitudes 

Love  waves 
percent  variation 
with  azimuth  (Icj) 

Rayleigh 

waves 

(vertical) 

Rayleigh  waves 
percent  variation 
with  azimuth  (la) 

Uniform 

0 

0 

136 

0 

Lithosphere 

5.5 

16 

137 

2 

Modified  lith. 

5.5 

18 

127 

2 

Lith.  (weakened) 

16.4 

17 

209 

3 

Mod.  lith.  (wknd) 

17.9 

19 

192 

3 

Note. — RMS  amplitudes  are  calculated  from  20  s  of  data  enclosing  the  waveforms  (67-87  s  for  Love  waves,  80-100 
s  for  Rayleigh  waves).  Waveforms  are  low-pass  filtered  at  2  Hz.  RMS  amplitudes  are  for  a  250  km  range  and  are 
scaled  by  le-7.  The  percent  variation  in  wave  amplitude  is  calculated  from  12  directions  (4  directions  in  each  of  3 
independent  calculations)  for  the  “lithosphere”  calculations  and  4  directions  from  the  “modified  lithosphere” 
calculations. 


^  iPun:  Iith1^  Direction:  W 


iiPun:  Iith1,  Direction:  W 


o  I 

T.  m  70  80  90  100 

Time 


^  tcPun :  lith  1 ,  Direction :  S 


Figure  43.  Surface  wave  displacements  in  4  directions  at  250  km  range  for  the  “lithosphere  1” 
calculation.  The  waveforms  are  low-pass  fdtered  at  2  Hz.  The  j^-axis  is  measured  in  units  of  le-6 
m  for  horizontal  displacements  and  le-7  m  for  vertical  displacements. 
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Figure  44.  Same  as  Figure  43,  but  for  the  “lithosphere  1  (weakened)”  calculation. 
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The  power  spectral  density  of  the  surface  waveforms  is  shown  in  Figure  45.  Similar  to  the  body 
waves,  there  are  nearly  negligible  differences  between  the  spectra  of  the  “lithosphere”  and 
“modified  lithosphere”  calculations.  On  the  other  hand,  the  inclusion  of  enhanced  material 
weakening  in  the  calculations  generates  a  more  noticeable  change.  In  the  case  of  body  waves  the 
spectra  have  maximal  power  near  ~4-8  Hz,  while  in  the  case  of  surface  waves  the  spectra  peak 
near  ~0.4-1.0  Hz. 
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Figure  45.  Normalized  power  spectral  density  of  the  surface  waves  for  displacements  in  the 
horizontal,  vertical,  and  radial  planes  (top  to  bottom).  The  right  panel  shows  the  power  spectra 
after  smoothing  in  frequency,  created  by  averaging  the  power  in  adjacent  frequency  bins  to  more 
clearly  show  the  spectral  differences  between  the  waveforms. 


The  difference  between  the  spectra  with  and  without  enhanced  weakening  is  presented  in  Figure 
46.  As  with  body  waves,  the  high  frequency  spectra  are  relatively  unchanged.  For  displacements 
along  all  three  directions,  material  weakening  increases  the  power  at  frequencies  below  ~0.5  Hz 
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and  between  ~0.8-1.5  Hz,  and  deereases  the  power  at  frequeneies  between  ~0.5-0.8  Hz.  Above 
1.5  Hz,  the  change  in  the  spectra  for  horizontal  displacements  relative  to  vertical  and  radial 
displacements  are  opposite  in  sign.  Horizontal  displacements  are  enhanced  from  1.5  Hz  to  3  Hz, 
and  are  reduced  from  ~3-7  Hz.  Vertical/radial  displacements  are  reduced  from  1.5  Hz  to  3  Hz, 
and  are  enhanced  from  ~3-7  Hz. 


Figure  46.  Difference  in  normalized  spectral  power  for  the  surface  waves  with  and  without 
enhanced  material  weakening. 

4.3.3  Heterogeneous  Source  Region:  Nonlinear  Deformation 

In  uniform  media,  symmetry  ensures  that  nonlinear  deformation  generated  by  the  explosion  is 
axisymmetric  along  any  horizontal  plane  (Figure  47).  In  this  case,  yielding  is  confined  to  a 
spherical  region  ~  150-200  m  surrounding  the  shot  point,  in  addition  to  a  broader  area  near  the 
surface  resulting  from  spall.  With  heterogeneous  material  properties,  the  axisymmetry  is  broken 
which  enhances  the  nonlinear  deformation  (Figure  48  and  Figure  49).  Rather  than  being  confined 
to  the  surface  and  a  spherical  volume  surround  the  explosion,  the  non-linear  deformation 
generated  by  an  explosion  in  heterogeneous  media  extends  through  ~10-100  m  filaments  from 
the  shot  region  to  the  surface.  In  some  cases,  heterogeneity  generates  channels  of  deformed 
material  emanating  up  to  400  m  from  the  source  region  (e.g..  Figure  48,  middle-right). 
Weakening  of  the  material  properties  has  a  marginal  (enhancing)  effect  on  the  extent  of 
nonlinear  deformation. 
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X(m) 


Figure  47.  Vertical  plane  through  the  shot  point  showing  nonlinear  damage  in  the  source  region 
comprised  of  uniform  material  properties.  Cells  that  have  yielded  are  shown  in  blue,  and  cells 
that  have  cracked  are  shown  in  red. 


X-Z  slices  Y-Z  slices 


Figure  48.  Top  and  middle  panels:  Vertical  plane  through  the  shot  point  showing  nonlinear 
deformation  in  the  source  region  for  the  “lithosphere  1”  and  “lithosphere  1  (weakened)” 
calculations,  respectively.  Bottom  panel:  Vertical  plane  through  the  shot  point  showing  the  shear 
modulus. 
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Figure  49.  Same  as  Figure  48,  but  for  the  “modified  lithosphere  1”  calculations. 

When  enhanced  weakening  is  included  in  the  calculations,  the  weakest  material  regions  exhibit 
greater  nonlinear  deformation.  Figure  50  shows  regions  of  enhanced  weakening  relative  to 
regions  of  nonlinear  deformation  for  two  of  the  calculations.  In  each  case,  the  area  of  nonlinear 
deformation  has  an  irregular  outline,  but  some  of  the  larger  protrusions  extend  into  the  weaker 
regions  where  they  abut  a  stronger  patch  of  material. 
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Enhanced  Weakening  (Lithosphere)  Enhanced  Weakening  (Modified  lith.} 


Figure  50.  Regions  of  enhaneed  weakening  (smoothed  with  a  Gaussian  kernel  of  20  m  half¬ 
width)  for  the  “lithosphere  1”  ealculation  (left)  and  the  “modified  lithosphere  1”  caleulation 
(right)  in  the  x-y  plane  at  shot  depth.  Overlaid  on  top  of  the  images  are  eontours  enelosing  the 
regions  of  nonlinear  deformation.  Eaeh  image  is  600  m  on  a  side.  Colors  span  0.35  to  0.65, 
indicating  the  fractional  strength  of  the  material  relative  to  the  non-weakened  calculations. 

5.  CONCLUSIONS 

In  this  project  we  are  trying  to  understand  the  effects  of  3D  structural  variations  on  seismic 
waves  from  underground  explosions.  One  of  the  main  objectives  is  to  understand  why  shear 
waves  generated  by  underground  explosions  are  ubiquitous.  We  have  performed  calculations 
with  variable  strength  distribution,  topography  and  3D  tectonic  strain  release,  and  modeled  the 
SHOAL  and  NPE  explosions.  We  have  run  a  calculation  of  the  NPE  with  and  without 
topography,  and  extended  the  calculation  to  400  km  using  the  representation  theorem.  The 
calculation  shows  that  the  effect  of  topography  can  be  modeled  and  produces  effects  consistent 
with  the  data.  Near  field  tangential  motion  is  caused  primarily  by  topographic  gradients.  The 
effects  of  near-source  topography  on  subsurface  ground  motion  is  primarily  due  to  off-axis 
propagation.  The  effects  of  near-source  topography  on  far- field  seismic  radiation  are  small. 

The  calculation  of  SHOAL  shows  that  tectonic  strain  release  consistent  with  the  local  stress  state 
can  generate  near  field  and  regional  signals,  including  long  period  surface  waves,  similar  to  those 
observed.  We  find  that  tectonic  release  causes  small  changes  to  the  far-field  P-wave  waveform, 
but  has  very  little  effect  on  far-field  P-wave  amplitudes.  It  does  generate  SH  waves  not  present  in 
the  calculation  without  tectonic  release.  This  implies  that  while  tectonic  release  may  significantly 
affect  Ms,  the  effect  on  mb  will  be  small.  The  calculated  effects  of  tectonic  release  are  similar  to 
those  observed,  but  not  quite  as  large.  This  is  probably  means  that  the  in-situ  material  which  is 
nearly  critically  balanced  against  tectonic  stresses  is  initially  weaker  than  the  material  model 
used  in  the  calculation  which  is  stronger  until  failure  occurs.  An  initially  weaker  material  would 
amplify  the  effects  that  we  observed  in  our  calculations.  It  is  very  likely  that  the  observed  near 
field  motion  and  Love  waves  and  Rayleigh  wave  variations  were  due  to  tectonic  strain  release. 

The  calculations  with  3D  heterogeneity  show  that  both  elastic  heterogeneity  and  strength 
heterogeneity  cause  generation  of  shear  waves,  in  particular  SH  motions,  however  the  effect  is 
much  stronger  when  strength  variation  is  included.  Strength  variability  leads  to  odd-shaped 
regions  of  nonlinear  deformation  that  enhance  shear  wave  generation. 
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